function hfile = glm_WriteAnalyzeBetas(hfile, tfolder)
% GLM::WriteAnalyzeBetas  - write beta images as Analyze
%
% FORMAT:       glm.WriteAnalyzeBetas(tfolder [, cv [, aformat]])
%
% Input fields:
%
%       tfolder     target folder
%
% No output fields.
%
% Note: RFX GLMs

% Version:  v0.8a
% Build:    9111921
% Date:     Nov-19 2009, 9:15 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
   ~ischar(tfolder) || ...
    isempty(tfolder) || ...
    exist(tfolder(:)', 'dir') ~= 7
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
tfolder = tfolder(:)';
if ~bc.ProjectTypeRFX
    error( ...
        'BVQXfile:InvalidObject', ...
        'Method only valid for RFX GLM objects.' ...
    );
end

% get subject and predictor names
subnames = regexprep(glm_Subjects(hfile), '[ \!\@\#\$\%\&\*\(\)\[\]]', '');
prednames = glm_SubjectPredictors(hfile);

% data handling depends on project type
imdims = size(bc.GLMData.RFXGlobalMap);
switch (bc.ProjectType)
    case {0}
        imsize = imdims;
        ires = 3;
        istart = 128 - round(1.5 * imdims);
    case {1}
        imsize = imdims;
        ires = bc.Resolution * ones(1, 3);
        istart = 128 - [bc.ZEnd, bc.XEnd, bc.YEnd] + ires;
    case {2}
        if imdims(1) > 50000
            ires = 2;
        else
            ires = 3;
        end
        
        % try to find good shape
        imsize = [2, 2, 2] * ceil((imdims(1) ^ (1 / 3)) / 2);
        istart = 128 - round(ires * 0.5 * imsize);
        imdims = imsize;
    otherwise
        error( ...
            'BVQXfile:InvalidObject', ...
            'ProjectType %d not supported.', ...
            bc.ProjectType ...
        );
end

% get Analyze framework right
ai = BVQXfile('new:hdr');
aibc = bvqxfile_getcont(ai.L);
aibc.FileMagic = 'ni1';
aibc.NIIFileType = 1;
aibc.ImgDim.Dim = [4, imdims, 1, 0, 0, 0];
aibc.ImgDim.DataType = 16;
aibc.ImgDim.BitsPerPixel = 32;
aibc.ImgDim.PixelSpacing = [0, ires, ires, ires, 1, 0, 0, 0];
aibc.DataHist.NIftI1.QFormCode = 2;
aibc.DataHist.NIftI1.SFormCode = 2;
aibc.DataHist.NIftI1.QuaternionB = 0;
aibc.DataHist.NIftI1.QuaternionC = 1;
aibc.DataHist.NIftI1.QuaternionD = 0;
aibc.DataHist.NIftI1.QuatOffsetX = istart(1);
aibc.DataHist.NIftI1.QuatOffsetY = istart(2);
aibc.DataHist.NIftI1.QuatOffsetZ = istart(3);
aibc.DataHist.NIftI1.AffineTransX = [ires(1), 0, 0, istart(1)];
aibc.DataHist.NIftI1.AffineTransY = [0, ires(2), 0, istart(2)];
aibc.DataHist.NIftI1.AffineTransZ = [0, 0, ires(3), istart(3)];
aibc.DataHist.NIftI1.NIftIMagic = 'ni1';

% preset VoxelData
if bc.ProjectType == 2
    aibc.VoxelData = single(zeros(prod(imsize), 1));
end

% iterate over subjects, predictors
for sc = 1:bc.NrOfSubjects
    for pc = 1:bc.NrOfSubjectPredictors

        % depends on Project Type
        switch (bc.ProjectType)
            case {0}
                aibc.VoxelData = bc.GLMData.Subject(sc).BetaMaps(:, end:-1:1, :, pc);
            case {1}
                aibc.VoxelData = permute( ...
                    bc.GLMData.Subject(sc).BetaMaps(end:-1:1, end:-1:1, end:-1:1, pc), ...
                    [3, 1, 2]);
            case {2}
                aibc.VoxelData(1:imdims(1)) = ...
                    bc.GLMData.Subjects(sc).BetaMaps(:, pc);
        end

        % set content and write as Analyze
        aibc.HdrKey.DBName = sprintf('BVQXGLM S%03d P%03d', sc, pc);
        aibc.DataHist.Description = ...
            sprintf('beta - sub %s, cond %s', subnames{sc}, ...
            prednames{pc});
        bvqxfile_setcont(ai.L, aibc);
        try
            aft_SaveAs(ai, sprintf('%s/%s_P%03d_beta.hdr', tfolder, subnames{sc}, pc));
            hdr_SaveVoxelData(ai);
            aibc = bvqxfile_getcont(ai.L);
        catch
            bvqxfile_clear(ai.L);
            rethrow(lasterror);
        end
    end
end

% clear object
bvqxfile_clear(ai.L);
